Properties of patchy colloidal particles close to a surface: a Monte Carlo and density 

functional study 
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We investigate the behavior of a patchy particle model close to a hard-wall via Monte Carlo 
simulation and density functional theory (DFT). Two DFT approaches, based on the homogeneous 
and inhomogeneous versions of Wertheim's first order perturbation theory for the association free 
energy are used. We evaluate, by simulation and theory, the equilibrium bulk phase diagram of the 
fluid and analyze the surface properties for two isochores, one of which is close to the liquid side of 
the gas-liquid coexistence curve. We find that the density profile near the wall crosses over from a 
typical high-temperature adsorption profile to a low-temperature desorption one, for the isochore 
close to coexistence. We relate this behavior to the properties of the bulk network liquid and find 
that the theoretical descriptions are reasonably accurate in this regime. At very low temperatures, 
however, an almost fully bonded network is formed, and the simulations reveal a second adsorption 
regime which is not captured by DFT. We trace this failure to the neglect of orientational correlations 
of the particles, which are found to exhibit surface induced orientational order in this regime. 



I. INTRODUCTION 

Colloidal particles with patterned surfaces - better 
known as patchy particles PHI) - have been studied ex- 
tensively in recent years owing to their ability to self- 
assemble in a rich number of cluster, gel, glassy and 
crystalline phases. Understanding how the surface pat- 
tern influences the self-assembly mechanism is crucial to 
a bottom-up strategy for designing new materials where 
the desired macroscopic behavior is encoded in the micro- 
scopic properties of the building-blocks Q. Patchy parti- 
cles represent a valuable model system for investigating 
and understanding the behavior of more complex con- 
stituents such as amphiphilic molecules, colloidal clays, 
proteins and DNA nano- assemblies New con- 

cepts, as equilibrium gel optimal network density [TlT | 
and empty liq uid have arisen from the study of the phase 
diagrams Il2l . [l3| of patchy particles with limited va- 
lence [l4|, Il5| . These models emphasize the role of the 
number of bonds between the particles in determining 
the equilibrium as well as the static and dynamic behav- 
ior of the system in and out of equilibrium . 

Homogeneous patchy particle fluids are described sat- 
isfactorily by Wertheim's first order thermodynamic per- 
turbation theory [l7l - f20| that provides an expression for 
the free energy per particle with n p patches, which are 
treated independently. Less known are the properties 
of patchy particle fluids in confined geometries. Under- 
standing the behavior of patchy particles close to surfaces 
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has direct impact on a number of different applications 
which require patchy particles to self-assemble in con- 
fined geometries. We recall for instance the templated 
self-assembly technique [22| where confined geometries 
are used to orient bulk structures or to induce the forma- 
tion of novel morphologies. Moreover recent studies have 
focused on effective forces between colloids generated by 
confined critical patchy particles, aiming to control col- 
loids stability [2l|. Finally such systems may be stud- 
ied with small-angle neutron scattering (23| and atomic 
force microscopy [24| and thus a quantitative microscopic 
description of patchy particles in confined geometries is 
highly desirable. 

The description of confined associating fluids has been 
addressed in the past. Density functional theory (DFT) 
based on a perturbation of the inhomogeneous hard- 
sphere fluid was used to describe the structure of inho- 
mogeneous associatin g fl uids. An early example is the 
work of Segura et al. [25| for particles with four patches 
close to a hard wall where the weighted density approxi- 
mation (WDA) of Tarazona [26|, [27j was combined with 
Wertheim's theory to obtain a DFT description of as- 
sociating fluids. In the original DFT formulation both 
homogeneous and inhomogeneous versions of Wertheim's 
first order perturbation theory were used to account for 
the association free energy of the particles and it was 
found that the homogeneous theory - where the law of 
mass action is identical to that of the bulk system with 
the density replaced by the averaged local density - yields 
satisfactory results over the whole range of parameters. 
By contrast, the inhomogeneous version of Wertheim's 
theory was found to overestimate badly the layering of 
the particles near the wall. The reasons for this dis- 
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crepancy are not fully understood but they appear to 
be related to the difficulty of implementing Wertheim's 
law of mass action in the inhomogeneous system. This 
approach was extended to mixtures of associating and 
neutral equi-sized hard- spheres [28[. 



A second approach for the same model - particles 
with four patches - was developed by Yu and Wu (29j . 
They employed a modified fundamental measure theory 
(FMT) [30| based on Rosenfeld's functional for inhomo- 
geneous hard-sphere fluids and an inhomogeneous ver- 
sion of Wertheim's free energy for association. Yu and 
Wu [lOl employed an inhomegeneous version of the law of 
mass action which was shown to give comparable results 
to the homogeneous theory of Segura et al. for interme- 
diate densities. Yu and Wu's theory, however, is more 
accurate at the highest densities and can be extended 
to binary mixtures of particles with different sizes. More 
importantly, it reveals that the inhomogeneous version of 
Wertheim's first order perturbation theory depends cru- 
cially on the generalization of the law of mass action for 
inhomogeneous patchy particle systems. 



Both theories were applied to systems where bonding 
is not fully developed, i.e. at temperatures that are not 
too low, and the relation between surface and bulk prop- 
erties has not been investigated. In this paper we fo- 
cus on the self-assembly of patchy particles with n p = 3 
at a planar hard wall, and extend our calculations to 
the region where a fully bonded optimal network devel- 
ops. We have simulated the system by fixing the density 
and scanning the temperature down to very low temper- 
atures. We compare the simulation results to the results 
of two density functional theories: a homogeneous WD A 
approach similar to that of Segura et al. and the inhomo- 
geneous FMT approach of Yu and Wu, which were shown 
to be equivalent for hard-spheres in the range of densi- 
ties considered here. For the three-patch particle model 
the theories give similar results at moderate to high tem- 
peratures, by contrast to the low-temperature regime, 
relevant to the formation of arrested states, where both 
theories break down. We trace this failure to the assump- 
tion of independent patches and stress the need for the 
development of density functionals that account for the 
directionality of the bonding, which plays a crucial role 
in the low- temperature regime. 




FIG. 1. Phase diagram of patchy particles with valence 
n p — 3. Symbols are results from MC simulations. Squares 
depict the gas-liquid coexistence line. Diamonds indicate clus- 
ter fluid states (pj, < 0.5), while triangles indicate percolating 
fluid states. Solid lines are the results from Wertheim's the- 
ory for the coexistence line and from Flory-Stockmayer theory 
for the percolation line. Dashed lines indicate the isochores 
investigated in this study. Patchy particles are modeled as 
hard spheres with 3 equidistant bonding sites (patches) on 
the particle equator (inset). 




FIG. 2. Comparison of the bonding probability pb from Monte 
Carlo simulations (symbols) and Wertheim's theory (Eq. [5] 
lines) . 



The manuscript is organized as follows: in Sec. II we 
describe the model and the simulation techniques. More- 
over, a brief comparison of the results of Wertheim's the- 
ory and the simulation results for the bulk fluid are pre- 
sented and discussed. In Sec. Ill we describe the two 
density functional approaches and in Sec. IV we discuss 
their validity and limitations, with emphasis on the low- 
temperature regime. 



II. BACKGROUND 
A. Model and simulation methods 

Patchy colloidal particles are modeled as hard spheres 
(HS) of diameter er with n p = 3 equidistant bonding sites 
on the particle equator (see the inset of Fig [I| . The in- 
teraction between two patches Vij belonging to particles 
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i and j is given by the Kern-Frenkel potential |3l| : 



B. Wertheim Theory 



Vu 



Vsw{\rij\)G(rij,h,rj), 



(1) 



where fij is the vector between the centers of particles i 
and j, and fi is the unit vector from the center of particle 
i to the center of one patch on its surface. Vsw is a square 
well potential, 



Vsw(\ni\) 



co if \fij \ < er, 
— e if er < If 
otherwise 



ij | < o 



(2) 



and G is the angular part: 



G(fij,fi, fj) 



1 if 



ij 

-f, 



Ti > COS 



rj > cos( 



(3) 



otherwise. 



The interaction energy between sites e sets the energy 
scale. The spatial range S and the angle 9 max control the 
volume available for bonding, Vb, which is: 

3 

Vb = ^~ K 1 + 5 l a ? -!][!- cos(6 max )} 2 . (4) 

We fix the parameters 5 = 0.119<7 and cos9 ma x = 0.895, 
fulfilling the single bond per patch condition assumed in 
Wertheim's first order perturbation theory. 

We perform Gibbs ensemble Monte Carlo (GEMC) (H 
simulations to locate the gas-liquid coexistence line and 
grand-canonical Monte Carlo (GCMC) simulations [33| 
to estimate the critical point. In the GEMC method, 
a MC step consists on average of 4000 roto-translation 
attempts, 400 particle swap attempts and one volume 
change. About 300 particles in a volume V = 2880er 3 
were simulated. In the GCMC method, we consider boxes 
with L = 6cr to L = 14a, and a MC step consists of 
500 roto-translation attempts and one insertion/deletion 
attempt. For the largest box, the number of particles 
fluctuates between zero and 800. We scan the chemi- 
cal potential /i and temperature T to locate the region 
where the system exhibits large fluctuations in the num- 
ber of particles N and in the energy E, which signal the 
presence of a critical point. The appropriate combina- 
tion of these fluctuations, at the critical point, follows 
the order parameter distribution of the Ising universality 
class (33.135). 

The surface properties are investigated by canonical 
ensemble Monte Carlo (MC) simulations, at different 
temperatures and densities, by introducing a surface in 
the middle of the simulation box and periodic boundary 
conditions. The surface is modeled by a planar hard wall 
acting on the particles through a hard-core repulsion, i.e., 
the interaction between the particles and the surface is 
purely entropic. The density profiles close to the wall 
are calculated at two bulk densities, at temperatures T 
above and below the percolation line. 



Within Wertheim's theory, the free energy for a homo- 
geneous system is the sum of two contributions. The free 
energy of the hard-sphere reference system and the bond- 
ing contribution Fb on d{p,T) which arises from consider- 
ing certain graphs in the Mayer expansion |3(|. For the 
present single-component model, the bonding contribu- 
tion can be expressed in terms of the bonding probability 
Pb (fraction of bonded sites) as 



/3F bond /N = n p ln(l 



, 1 
Pb) + -^n p p b , 



(5) 



where j3 = 1/kT with k the Boltzmann constant. Assum- 
ing that all sites have the same probability of bonding, 
Wertheim's theory predicts that pb is determined by the 
law of mass action: 



Pb 



{i-Pb? 



pn p A, 



(6) 



where A is the equilibrium constant of the "reaction" 
(bonding) between two patches and p is the density. To 
evaluate A for the present model, we assume that the ra- 
dial distribution function of the reference HS fluid gnsi 7 ") 
is approximated by (3^ | : 



9Hs(r) = (A + Ax) + Atir/a - 1), 



where 



A = 
Ai = 



As a result, 
A = 4ttx 2 



I-O.577 4.577(1 + ??) 

-4.577(1 + 77) 

(1 — v) 3 ■ 

(1 + 5) 3 - 1 



. (l + ,5) 4 -l 
A + f A x 



(7) 

(8) 
(9) 

(10) 



x [exp(/3er) - 1] . 

is the packing fraction and x 



TV pa 



0.5(1 



where 77 — 

cos(9 rnax )) is the fraction of surface covered by the patch. 
Note that the linear approximation of Eq. [7| is highly 
accurate in the relevant r range, i.e. within the well of 
the square-well potential. 

The phase diagram is then calculated straightfor- 
wardly. It includes a gas-liquid first order phase tran- 
sition that ends at a critical point at low T and p. The 
coexisting homogeneous phases have different densities 
and fractions of unbonded sites. The percolation line 
is calculated using the Flory-Stockmayer (FS) theory of 
polymerization |38| - |40| | which gives for the percolation 
threshold, 
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Pb = 



(11) 
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C. Bulk behavior 

Fig. Q] shows the phase diagram obtained via GCMC 
and GEMC simulations and compares it with the theo- 
retical results from Wertheim's theory. The figure also 
shows the FS percolation line. The high density phase is 
always percolated, in the sense that there is a non-zero 
probability of finding an infinite cluster that contains al- 
most all the particles. The percolation line (i. e. the line 
that separates percolated from non-percolated states) in- 
tercepts the binodal on the low-density phase, near the 
critical point. As shown in Fig. [TJ the percolation line 
starts on the left of the critical point and the temper- 
ature increases monotonically with the density. In Fig. 
[21 we plot the bonding probability as a function of tem- 
perature for different densities. The agreement between 
Wertheim's theory and Monte Carlo simulations is quite 
good both for the phase diagram and for the bonding 
probability. 

III. DENSITY FUNCTIONAL THEORY 

As usual in density functional theory, we split the 
Helmholtz free-energy functional into the ideal and ex- 
cess parts: 

F[p(r)} = F ld W)] + JWp(r)], (12) 

where p(r) is the number density. Here and in what 
follows, we assume that the single-particle distribution 
function p(r) depends on the spatial but not on the ori- 
entational coordinates. Note that the arrangement of 
patches on the equator breaks the rotational symmetry 
of the particles, and therefore states with a preferred ori- 
entation of the particles cannot be ruled out. In fact, 
this is a crude approximation that will break down at 
low temperatures, as we will discuss later. 
The ideal part is given by: 

f3T ld [p{r)] = J d 3 rp(r) [ln(p(r)V) - 1] , (13) 

where V is the thermal volume. The integral is over the 
volume V. 

T exc contains the excluded volume interactions be- 
tween HS and the bonding free energy due to bond for- 
mation between the particles. We have used two different 
approximations for J- exc : a modified version (4l| of the 
local weighted density approximation (WDA) introduced 
by Segura et al. (25| and the fundamental-measure the- 
ory (FMT) for associating fluids developed by Yu and 
Wu [H]. 

In the WDA method, the HS and the association term 
in the free-energy per particle are evaluted for a homo- 
geneous system at the same weighted density. The lat- 
ter is calculated for a fluid of HS |4l|. The Carnahan- 



Starling (36| approximation is employed for the HS con- 
tribution, while the association free energy is given by 
Wertheim's first-order perturbation theory for a homo- 
geneous system(Eq. [5j. 

In the FMT proposed by Yu and Wu [H|, Wertheim's 
free-energy functional for the inhomogeneous system is 
used and thus the two contributions are treated sepa- 
rately. The HS reference fluid is described by the Rosen- 
feld FMT approach while the association contribution 
is given by an appropriate inhomogeneous Wertheim's 
term. A detailed description of both methods is found in 
appendices [A] and [B] respectively. 



IV. RESULTS 

We focus on the surface properties of the fluid at two 
bulk densities (pb<7 3 = 0.70 and p&o -3 = 0.40) in contact 
with a neutral hard- wall, and consider different tempera- 
tures (see Fig. [TJ . The main results are reported in Fig. [3] 
and Fig. [H 

The density profiles calculated from MC simulations 
for pbcr 3 = 0.70, are plotted together with FMT and 
WDA theoretical results in Fig. [3] The limiting case 
of HS is also included for reference (panel (a) of Fig. [3J . 
Panels (b)-(d) of Fig.[3Jshow the results for three different 
reduced temperatures T above and below the percolation 
line. The main difference between the associating fluid 
(finite temperature) and the reference HS fluid (T — > oo) 
is the reduction in the adsorption of particles near the 
wall. This effect can be understood as follows: at very 
high temperatures almost all the particles are fully un- 
bonded, and there is a strong adsorption of particles near 
the wall due to entropic reasons (gain in configurational 
entropy). As we decrease the temperature, the level of 
association between particles increases. The probability 
of bonding at distances of order a from the wall is lower 
than at larger distances (the wall is neutral). The con- 
tribution to the bonding free energy of the particles near 
the wall is lower than the contribution due to the other 
particles, resulting in a reduction of the adsorption of 
patchy particles at the wall when compared to the ad- 
sorption of hard-spheres (see the contact density, i. e. 
the value of the density at z/a = 0.5, in panel (f) of Fig. 

ej. 

The WDA provides a fairly good description of the 
density profiles down to kT/e = 0.15, i.e. well inside 
the percolation region (see Fig. [TJ. At this temperature 
WDA seems slightly more accurate than FMT, providing 
a good estimate of the contact density. The agreement is 
lost at kT/e = 0.10 and below (Fig.[3](d) and (e)) where 
both theories underestimate the contact density and do 
not account for the strong fluid layering on increasing 
the distance from the wall. At this density, we find that 
patchy particles are adsorbed on the wall at all T (al- 
though the adsorption decreases as the temperature de- 
creases). No changes with further cooling are expected at 
lower temperatures since pi,(kT/e — 0.08) ~ 0.98, mean- 
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FIG. 3. Density profiles as a function of the distance from the wall at different temperatures for pba :i — 0.70: (a) kT/e — > oo, 
(b) kT/e = 0.30, (c) kT/e = 0.15, (d) kT/e = 0.10, and (e) kT/e = 0.08. The full circles are Monte Carlo simulation results. 
The lines are the results from density functional theory: FMT (red solid line) and WDA (dashed blue line) . (f ) Contact density 
as a function of temperature. 



ing that the liquid has reached an almost fully bonded 
configuration and thus the structural properties become 
essentially T-independent 0, 53] ■ 

In Fig. [4l panels (a)-(g), we plot the density profiles 
for phcr 3 = 0.40. Panel (a) illustrates the HS fluid corre- 
sponding to T — > oo. At kT/e = 0.3 (b) and kT/e = 0.2 
(c) both, FMT and WDA, are in reasonable agreement 
with the MC simulation results. FMT yields slightly bet- 
ter results than WDA for kT/e = 0.2. Again, the two 
density profiles indicate an adsorption of particles near 
the wall. The adsorption, due to the excluded volume, 
is small compared to the adsorption of HS (Fig. H (a)). 
Fig. g] (d) illustrates the system at kT/e = 0.15 (slightly 
below the percolation threshold). The MC and FMT 
density profiles are almost uniform. Patchy particles are 
slightly desorbed from the wall, indicating the cancel- 
lation between the effects due to excluded volume and 
association. At this temperature WDA predicts a des- 
orption of particles very close to the wall, but it also pre- 
dicts the formation of a more or less well-defined layer of 
particles at approximately z/a = 1 from the wall. Pan- 
els (e)-(g) illustrate the low temperatures: kT/e = 0.12 
(e), kT/e = 0.10 (f), and kT/e = 0.08 (g). At these T 
the association of particles is very high (see the fraction 
of unbonded sites in Fig. [5]), and the energy of bond- 
ing dominates the behaviour of the system. As a result, 



FMT predicts a strong desorption of particles from the 
wall. However, the simulation shows the opposite behav- 
ior; there is a well-defined layer of particles near the wall. 
The layer grows and approaches the wall as the temper- 
ature decreases. This desorption-adsorption crossover is 
also reflected in the contact density depicted in panel (h) 
of Fig. |U WDA predicts the presence of a layer of parti- 
cles at approximately z « a for low T. Nevertheless, the 
peak heights do not vary significantly with T, by contrast 
to the simulation results. 

The overall T dependence of the contact density at 
Pbcr 3 = 0.40 is significantly different from the system at 
Pbcr 3 = 0.70. Indeed, for p^a 3 = 0.40, the contact den- 
sity changes continuously from the HS adsorption limit 
to the low T desorption, while in the p^cr 3 = 0.70 system 
only adsorption is present. To rationalize this behavior 
we recall the thermodynamics of the bulk system, and in 
particular the gas-liquid coexistence. In limited valence 
systems, the liquid side of the gas-liquid coexistence is al- 
most vertical in the T — p plane. The associated density 
provides a quantification of the so-called optimal network 
density, i.e. the density at which particles in the liquid 
(actually a gel at low T) are not stressed. At the same 
time, the small value of the density of the gas-phase in- 
dicates that the coexisting pressure is rather small. The 
contact density, at a hard-wall, is a direct measure of the 



(i 




0.15 

kT/e 



FIG. 4. Density profiles as a function of the distance from the wall at different temperatures for p6<7 3 = 0.40: (a) T — s> oo, (b) 
kT/e = 0.30, (c) kT/e = 0.20, (d) kT/e = 0.15, (e) kT/e = 0.12, (f) kT/e = 0.10, and (g) kT/e = 0.08. The full circles are 
Monte Carlo simulation results. Lines are density functional theory results: FMT (red solid line) and WDA (dashed blue line), 
(h) Contact density as a function of temperature. 
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pressure in the system. Hence, close to the coexisting liq- 
uid branch, the contact density decreases as T decreases, 
resulting in the reported desorption. On increasing p at 
constant low T, the formation of an extended network 
of bonds causes the increase of stresses in the system 
and the pressure increases. This results in a significantly 
large value of the contact density, driving the adsorption 
phenomenon. 

To pin down the origin of the discrepancies between 
the DFT results and MC simulations which build up on 
cooling, we plot in Fig. [5] (a) the uniaxial order parameter 
profile Sn{z) for the system at pi,<T 3 = 0.40. This is de- 
fined as Sn(z) = J dCt h(Ct, z) P2(cos(#)) where h(Cl, z) is 
the orientational distribution function at distance z from 
the wall, P2(cos(0)) is the second Legendre polynomial 
and 6 is the angle between the unit vector u normal to the 
wall and the unit vector p normal to the plane containing 
the patches. Sn(z) provides information on the orienta- 
tion of particles as a function of the distance z from the 
wall. As shown in Fig. [SJ Sn(z) grows significantly close 
to the wall on cooling. The value of Sn(z) signals a 
cross-over from an isotropic fluid to a nematic-like phase 
near the wall. To highlight the particles orientation close 
to the wall, we show in Fig. \5\b) the probability distri- 
butions P(9,z), confirming that at low T the particles 
are oriented near the wall, with the plane containing the 
patches parallel to the wall. This geometry maximizes 
the bonding probability by moving the patches away from 
the neutral wall. The average orientational angle close to 
the wall decreases continuously with T and it should ap- 
proaches 9 = for perfect order. We note indeed that 
the finite bonding volume allows for a flexibility in the 
orientation of the particles, contributing to a small but 
non-zero average angle even in highly bonded conditions 
close to the wall. At high T the particles are randomly 
oriented (see Fig. [5] (a) and (b)) and P(9, z) ~ sin(0). Fi- 
nally we note that far from the wall (Fig. 0(c)) particles 
undergo a continuous change from a state with preferred 
orientation to an isotropic one. 

We stress that both DFT approaches neglect the ori- 
entational order that develops in the system at low T, 
since the free-energy depends on the number density but 
not on the orientation of the particles. 

Finally we comment on the difference between the re- 
sults of the two DFT approaches at low T for p&er 3 = 0.40. 
Even though the agreement with the numerical results is 
poor in both cases, WDA describes qualitatively the lay- 
ering of the particles near the wall at low T while this 
feature is absent in the FMT results. Such layering is re- 
lated - in part - to the orientation of the particles (which 
is neglected in both DFTs) in order to minimize the frac- 
tion of unbonded sites X(z) (i.e. f —pb) of particles near 
the wall. Fig.[5]shows X(z) along z as predicted by FMT 
(solid lines) and WDA (dashed lines). Symbols are the 
results from MC simulations. At low T, X(z) obtained 
from the homogeneous Wertheim theory (WDA) exhibits 
a trend similar to the MC results, while X(z) calculated 
using the inhomogeneous Wertheim theory (FMT) in- 
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FIG. 5. (a) Orientational order parameter Sn(z) near the wall 
extracted from MC canonical simulation of patchy particles 
for different T and pi,a 3 = 0.40. (b) Probability distributions 
P(9, z) as a function of 9 for z — 0.75a at the same T and p as 
in panel (a), (c) Probability distributions P(9, z) as a function 
of the distance from the wall z evaluated at kT/e = 0.1 and 
pa 3 = 0.4. 9 is evaluated at different distances from the wall 
by dividing the simulation box into "slices" of size L x L x Az 
where Az = 0.5a. Notice that when orientational isotropy is 
restored (high T), P(8) ~ sin(#). Lines are guides to the eye. 
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FIG. 6. Fraction of unbonded sites as a function of the nor- 
mal distance from the wall at different Ts and pb<J 3 = 0.40. 
Symbols are results from MC simulation. Solid lines are re- 
sults from the FMT-inhomogeneous Wertheim theory (eq. 
IB20) . Dashed lines are obtained from the WDA-homogeneous 
Wertheim theory (eq. [5]). 



the contact density at the hard wall is a measure of this 
pressure, at small pb desorption must occur. On increas- 
ing T the pressure increases and for T above the perco- 
lation threshold the particles are adsorbed at the wall. 
Such crossover is not observed at the higher pb where the 
density at the wall is always larger than the bulk density. 
Indeed, even at low T, the price to pay for the formation 
of a distorted network leads to an increase of the pressure 
and hence to a large contact density, even when pb — > 1 . 

Not surprisingly, the description using density func- 
tional theory is consistent with simulations and with the 
results reported in Ref. [HJ and [2!| at high and interme- 
diate T, but it fails at low T. We have traced this failure 
to the inadequacy of describing the orientational degrees 
of freedom of the particles. Indeed, close to the wall, the 
particles are oriented in such a way that the plane con- 
taining the bonding sites is almost parallel to the wall. 
It appears that the orientational degrees of freedom and 
an appropriate coupling to the density profile need to be 
taken into account in future work in order to describe the 
structure of associating fluids (and gels) close to a hard 
wall at low temperatures. 
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creases sharply near the wall. The WDA X(z) follows 
the density profile, increasing over the bulk value when 
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V. CONCLUDING REMARKS 

In this work we have studied the properties of a fluid 
of patchy particles with n p = 3 patches. We have located 
the gas-liquid coexistence line, the critical point and the 
percolation line, providing a full characterization of the 
thermodynamics and structure of the fluid phases. We 
have also shown that, for this patchy model, Wertheim's 
theory for homogeneous fluids is accurate. We have then 
studied the surface properties of this fluid in contact with 
a neutral hard-wall using MC simulation and two DFT 
approaches. 

We have investigated the behavior of the system at 
two different densities, one close to the liquid branch of 
the coexistence curve and one about 70 per cent higher, 
for several T, covering the structural change of the sys- 
tem from a monomer solution to an almost fully bonded 
network state. We have shown that at low pb, the wall 
adsorbs particles at high T which desorb as T is lowered. 
The physical mechanism responsible for the adsorption- 
desorption cross-over is understood in terms of the prox- 
imity to the gas-liquid coexistence curve. Indeed, at low 
T, close to the liquid branch, the liquid coexists with a 
gas at very low density and the pressure is small. Since 
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Appendix A: Weighted density functional 

The weighted density functional approach is based on 
the idea that the free energy of the inhomogeneous sys- 
tem, characterized by the single-particle density p(r), 
may be written in terms of the free energy of a homo- 
geneous fluid with an effective density that is evaluated 
through an appropriate averagin g pr ocedure. The ap- 
proach was proposed by Tarazo nal27i a nd was used and 
modified by several authors (30l. l43l. |44| . 

The excess free energy functional of the inhomogeneous 
fluid is written as 

Fe*c\p] = J drp^ftp;^, (Al) 

where f[p; r\ is the local excess free energy functional per 
particle. f[p; r\ can be written as a function of the den- 
sity p, a functional of the single-particle density, which 
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satisfies 



The grand potential functional 



/[p;f] = /(p(r)), 



(A2) 



where f(p(r)) is the local free-energy density of the ho- 
mogeneous system and p(r) is the weighted density de- 
fined by: 



p{r) = I dr 1 'p(r*)ui(r — ? '; p(r)). 



(A3) 



In Eq. IA3I p{r) is the difference between the density 
and its weighted counterpart and uj is the weight function 
satisfying the constraint 



dr uj(f — r , p(f)) = 1. 



(A4) 



The weight function, which encodes the non-local charac- 
ter of the functional, is related via non linear differential 
equations to the direct correlation function of the inho- 
mogeneous fluid. An approximation to the weight func- 
tion is obtained by requiring that the second functional 
derivative of the excess free-energy 



c (2, (ri,r 2 ;p) = - 



!35 2 T ex [p] 
5p{r 1 )5p{r 2 ) 



(A5) 



gives an accurate description of the correlations of the 
homogeneous fluid. 

Different WDA functionals result from different ap- 
proximations for the weight function. We use a modified 
version of the WDA proposed by Tarazona, developed by 
Kim and coworkers [4l| . In this approximation 



p{r) = / dr'p(r')uj(r-f',p b ), 



(A6) 



where p b is the bulk density. Following Tarazona, w(r,p) 
is expanded in powers of p in order to reduce the com- 
putational effort: 



ui(r, p) = u (r) + uJi(r)p + u 2 (r)p 2 . (A7) 
The same is done for p(r): 

P(r) = Po (r) + Pi (r)p(r) + 92 (f)p(f) 2 , (A8) 

where 



with 



Pi(r) 



p(r) = p a (r) + Pl (r)p b + p 2 {r)p 2 , (A9) 



df l p(r f )Lu l (f-r) i = 0,1,2. (A10) 



pn\p] = f3F[p] + /3J dr(V ext (f) - p)p{r) (All) 
is minimized with respect to variations of p(r): 



f35n[p] _ PSFjp] 
5p(f) Sp(r) 



-P0*-v ext (f)) = o> 



(A12) 



which yields the equation for the equilibrium density pro- 
file, 



PSF[p] 
8p{r) 



f3p-/3V ext (r). 



(A13) 



The first derivative of the excess free energy functional 
is the single-particle direct correlation function c^'{f) 
and putting together Eqs. [T31 IAU IA2I and IA3I we obtain 



Pfi - pV ext {r) = lnp(r) - c«[r;p], 



with 



c^[r;p\ = -f3f(p(r)) - (3 / d^p^)f'{p(f)) 



and 



(A14) 



Sp(r) ' 
(A15) 



cj(r- r',p(r')) 



5p(r) 

+ u(f-f,pb) J dr"p{f 4 ')oj'{f' -r",p(r*)). 

(A16) 

The chemical potential pi is evaluated from the homo- 
geneous version of Eq. IA15I 



j3p = lnp fc - c {1 \p b ) 

= lnp fc - 0f(pb) - (3p b f'(p b ). 



(A17) 



Finally, combining Eqs. lA14l and lA17l we obtain the den- 
sity profile: 

p(f) = p b exp[-^^t (r) + c« [r; p] ~ c« (p 6 )] . (A18) 

The excess free energy is the sum of the free energ y of 
HS given by the Carnahan-Starling approximation [36fl 
and the bonding contribution given by Eq. [SJ 

At a planar hard- wall, the external field is V ext = 
for z > and infinite otherwise, and the density profile 
depends only on the distance z from the wall: 
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p(z) = p 6 exp[c«(z; [p]) - c«(p b )] z > a/2, 

= z < a/2. (A19) 

The weight functions u>i(r) for this system, defined in 
Ref. Ell, are: 



Z < a, 



47TCT 3 ' 

otherwise, 



(A20) 



oj-l(z) = 0.21cr 2 - 1.49z 2 - 0.18( : 



1.36 , z < a 

a 



-O.llo- 2 + 2.9z 2 + 0.29(— Y - l.81a\z\ 

s 

z 3 

-1.6— a < z <2a, 
a 

otherwise, (A21) 



and 



IOtt^ 



a(a A - 12(crz) 2 - 5z 4 + 16a|^| 3 ), \z\ < a, 



576 

otherwise. 



(A22) 



Appendix B: Fundamental measure density 
functional 



PJhs = -no ln(l - n 3 ) + 

(n 2 ) 3 - 3n 2 n v2 ■ n v2 
24tt(1 - n 3 ) 2 



nin 2 - n v i ■ n v2 



1 - n 3 



(B2) 



where we have dropped the spatial dependence of the 
weighted densities for convenience. n a (r) are convolu- 
tions of the density with the weight functions w a (r), 
which are related to geometrical properties of the par- 
ticles: 



n a {r) = w a (f) *p(r), 



(B3) 



where * denotes the three-dimensional convolution h(r) * 
g(r) = J d 3 xh(x)g(x — r). For HS the weight functions 
are: 



w 3 (r) = e(R s -\r\), 
Mr) = HRs - H), 
wi{r) = W2(r)/{4nRs), 
w (r) = w 2 (r)/(4nR 2 s ), 

T 

w v2 {f) = w 2 (r) — , 

In 

w vl (f) = w v2 (r)/(4:irR s ). 



(B4) 
(B5) 
(B6) 
(B7) 

(B8) 

(B9) 



<$(•) is the Dirac-delta distribution and 9(-) is the Heav- 
iside step function. R$ = a/2 is the sphere radius. In 
planar geometry the one particle distribution function 
depends only on the normal distance from the wall, z. 
The weight functions are obtained integrating over the 
lateral coordinates, 



A second approach, uses a geometry-based density- 
functional to describe the excluded volume between 
spheres as well as the inhomogeneous associating free- 
energy. It was proposed by Yu and Wu (2j| to describe 
inhomogeneous mixtures of hard-spheres (HS) with an 
arbitrary set of interaction sites. In the following we 
consider a single-component fluid of HS with n p identical 
patches, in which case, the excess free-energy reads: 

P T exc \p{r)\ = fd 3 r {f3f HS [n a (r)] + f3f bond [n a (r )] } , 



(Bl) 

where fas is the reduced excess free-energy of the fluid 
of HS, and f bo nd is the free-energy arising from the asso- 
ciation of the particles. Both quantities depend on a set 
of weighted densities {n a (r}} (see later). 



w a {z) = J dx J dyw a {r). 
The resulting weight functions are 

w 3 (z) = ir(R 2 s -z 2 )e(R s ~\z\), 
w 2 (z) = 2nR s e(Rs-\z\), 
wi(z) = w 2 (z)/(4irRs), 
W (z) = w 2 (4irR 2 s ), 
w v2 = 2irzO(R s - \z\)z, 
w v \ = w v2 /{AttR s ), 

with z the unit vector normal to the wall. 



(BIO) 



(Bll) 
(B12) 
(B13) 
(B14) 
(B15) 
(B16) 



a. Hard sphere fluid 

The excluded volume interactions between HS are de- 
scribed by the Rosenfeld functional [3Cj: 



Wertheim's inhomogeneous free-energy 



The bulk free ener gy o f a fluid of particles with n p 
identical sites (eq. [S]) |17H20l | can be rewritten in terms 
of the bulk fraction of unbonded sites X b = 1 — p b as 
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Pfbond = n p p b In X, 



2 



(B17) 



with pi, the bulk density, and Xb the bulk fraction of 
unbonded sites, related to the thermodynamic variables 
through the law of mass action: 



X b = (1 + n pPb X b Ay 



(B18) 



Yu and Wu [29j generalized the bulk free energy, Eq. 
(|B17p , to inhomogeneous systems by including a new fac- 
tor £ = 1 — n v in v ijr^ that incorporates the vectorial 
weight densities into the associating part of the free en- 
ergy. The bonding free energy for inhomogeneous fluids 
reads: 

X(f) V 
2 2, 

(B19) 

where X{r) is the fraction of unbonded sites at position 
r given by the modified law of mass action: 



f3fbond[n a (r)] = n p n {r)C,(r) lnX(r) 



X{r) = (1 + n p n (r)ar)X(r)A(r)Y 



(B20) 



The interaction between two sites determines A. For 
the Kernel-Frenkel potential (where orientational and 
translational degrees of freedom are decoupled) A is given 



in terms of gHsi 1 ^)- the pair correlation function of the 
reference HS fluid, and Jm, the Mayer function: 



A(r) = / d 3 rg HS (r)f M , 



(B21) 



where /m — cxp(/3e) — 1, and the integral is over the 
bonding volume Vb- Following Yu and Wu (29| we use a 
modified contact value for the pair correlation function: 



gHs(r) 



1 



an 2 ( 



<x 2 (n 2 ) 2 C 



l-n 3 4(1 -n 3 ) 72(1 - n 3 f 



(B22) 



Assuming that the pair correlation function is constant 
over the bonding volume, we approximate Eq. (|B21[) by: 



A(f) = v b gHs(r)fM- 
Finally, we minimize the grand potential, 

Q\p] = F[p] -pi d 3 rp(r), 



(B23) 



(B24) 



to obtain the equilibrium density profiles, p is the chem- 
ical potential and the integrals are computed using the 
trapezoidal rule and a step size Az = O.Olcr. We use a 
standard conjugated-gradients method to minimize Q. 
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